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ABSTRACT 

We model the polarized thermal dust emission from protostellar cores that are assembled by 
£C) , super-sonic turbulent flows in molecular clouds. 

Self-gravitating cores are selected from a three dimensional simulation of super-sonic and 
super-Alfvenic magneto-hydrodynamic (MHD) turbulence. The polarization is computed 
in two ways. In model A it is assumed that dust properties and grain alignment efficiency 
are uniform; in model B it is assumed that grains are not aligned at visual extinction larger 
than Ayo = 3 mag, consistent with theoretical expectations for grain alignment mechanisms. 
Instead of using a specific set of grain properties, we adopt a maximum degree of polarization 
Pmax = 15 %. Results are therefore sensitive mainly to the topology of the magnetic field 
(model A) and to the gas distribution that determines the distribution of Ay (model B). 
Furthermore, the radiative transfer in the MHD model is solved with a non-LTE Monte Carlo 
method, to compute spectral maps of the J=l-0 transition of CS. The CS spectral maps are 
used to estimate the turbulent velocity, as in the observations. 

The main results of this work are: i) Values of P between 1 and 10% (up to almost P m ax) 
■ are typical, despite the super-Alfvenic nature of the turbulence; ii) A steep decrease of P with 

increasing values of the sub-mm dust continuum intensity I is always found in self-gravitating 
cores selected from the MHD simulations, if grains are not aligned above a certain value of visual 
extinction Ay.o (model B); iii) The same behavior is hard to reproduce if grains are aligned 
independently of Ay (model A); iv) The Chandrasekhar-Fermi formula, corrected by a factor 
/ w 0.4, provides an approximate estimate of the average magnetic field strength in the cores. 

Sub-mm dust continuum polarization maps of quiescent protostellar cores and Bok globules 
have recently been obtained. They always show a decrease in P with increasing value of / 
consistent with the predictions of our model B. We therefore conclude that sub-mm polarization 
maps of quiescent cores do not map the magnetic field inside the cores at visual extinction larger 
than Avfi ~ 3 mag. The use of such maps to constrain models of protostellar core formation 
and evolution is questionable. This conclusion suggests that there is no inconsistency between 
the results from optical and near-IR polarized absorption of background stars, and the observed 
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polarization of sub-mm dust continuum from quiescent cores. In both cases, grains at large 
visual extinction appear to be virtually unaligned. 

Subject headings: turbulence - ISM: kinematics and dynamics - radio astronomy: interstellar: 
dust continuum, polarization and lines 

1. Introduction 

The topology of the magnetic field inside and around protostcllar cores can be predicted with models 
for their formation and evolution. Models of sub-critical cores contracting under the effect of gravity and 
ambipolar drift, for example, predict a rather uniform magnetic field roughly perpendicular to the core 
major axis. Such models can be tested if the magnetic field topology is constrained observationally. 

Polarization maps of sub-mm thermal dust emission have recently been obtained for a number of 
protostellar cores and Bok globules (Minchin, Bonifacio & Murray 1996; Glenn, Walker & Young 1999; 
Greaves et al. 1999; Coppin et al. 2000; Davis et al. 2000; Matthews & Wilson 2000; Vallee, Bastien 
& Greaves 2000; Ward-Thompson et al. 2000; Henning et al. 2001). Results are still sparse, but are 
readily used to try to test theoretical models of core formation and evolution (Glenn, Walker & Young 
1999; Coppin et al. 2000; Davis et al 2000; Ward-Thompson et al. 2000). Since the interpretation of the 
observations is hardly unique, it is useful to compute polarization maps of the theoretical models, where 
the full three-dimensional information is available, and to infer their observable properties. Theoretical 
polarization maps have been previously computed by Wardle & Konigl (1990) for the molecular disk at the 
Galactic center and by Fiege & Pudritz (2000) for molecular cloud filaments. 

In this work we assume that protostellar cores are assembled by super-sonic turbulent flows in molecular 
clouds. We have shown in a previous work that the formation of cores by turbulent shocks provides an 
excellent interpretation of the relation between integrated intensity and rms velocity in molecular clouds 
(Padoan et al. 2001). Alternative interpretations of this newly discovered property of molecular clouds have 
not been proposed so far. Furthermore, super-sonic turbulence is ubiquitously observed in molecular clouds, 
and therefore models of core formation and evolution must take into account this turbulent environment 
self-consistently. We do so by selecting self-gravitating cores from a numerical simulation of super-sonic, 
super- Alfvenic, self-gravitating magneto-hydrodynamic (MHD) turbulence. The simulation is intended to 
represent a molecular cloud region with size L — 6.25 pc, rms turbulent velocity a v = 3.0 km/s and average 
gas density (n) = 320 cm~ 3 . These are typical values for molecular clouds according to Larson's relations 
(Larson 1981). The three selected cores have a mass of « 60 M Q and a maximum gas density ~ 10 5 cm~ 3 . 

The polarized thermal dust emission is computed in two ways. In model A it is assumed that dust 
properties and grain alignment efficiency are uniform; in model B it is assumed that grains are not aligned 
at visual extinction larger than Ay$ — 3 mag, in line with theoretical expectations for grain alignment 
mechanisms (Lazarian, Goodman & Meyers 1997). Instead of using a specific set of grain properties, we 
adopt a maximum degree of polarization P max = 15 %. Results are therefore sensitive mainly to the 
topology of the magnetic field (model A) and to the gas distribution that determines the distribution of Ay 
(model B). A Monte Carlo method is used to compute the distribution of Ay in the MHD model. 

We find that despite the relatively weak magnetic field averaged over the whole computational box, 
(B) = 3.3 ^G and (B 2 ) 1 / 2 = 7.1 /iG, the degree of polarization is typically between 1 and 10% (as 
observed), spanning a range of values almost up to P ma x- In protostellar cores (Davis et al. 2000) and 
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Bok globules (Hcnning ct al. 2000), P is always found to decrease with increasing sub-mm dust continuum 
intensity /. This is consistent with our MHD model only if grains are not aligned above a certain value 
of the visual extinction, Ay (model B), in agreement with previous results from optical and near-IR 
polarized absorption of background stars (Goodman et al. 1995; Gerakines et al. 1995; Arce et al. 1998). 

We also solve the radiative transfer in the MHD model with a non-LTE Monte Carlo method, to 
compute spectral maps of the J=l-0 transition of CS. We use these maps to estimate the turbulent 
velocity dispersion in the same way as in the observations. This turbulent velocity dispersion is used, 
together with the dispersion in the polarization angle, to infer the magnetic field strength using the 
Chandrasekhar-Fermi formula (Chandrasckhar & Fermi 1953). It is found that the volume-averaged three 
dimensional magnetic field strength can be estimated, despite the large dispersion in the polarization angle, 
if the Chandrasekhar-Fermi formula is used with a coefficient / 0.4. 

The polarization properties of turbulent molecular clouds depend on the average magnetic field 
strength, relative to the kinetic energy of turbulence, as measured by the rms Alfvenic Mach number, A^a, 
that is the ratio of the rms flow velocity and the Alfven velocity, (Clemens et al. 2000; Ostriker, Stone 
& Gammie 2001). A complete statistical study of the polarization properties of turbulent clouds with 
different values of M a will be presented separately, and will be more useful in the future, when a larger 
number of polarization maps of sub-mm and far-IR dust emission will be available (especially if future 
satellite missions of the type of M4 (Clemens et al. 2000) will be dedicated to this purpose). In the present 
work we focus on the properties of a few cores, to provide insight in the interpretation of recent sub-mm 
observations. 

The MHD models are presented in the next section, while in § 3 we describe the method for computing 
the polarized dust emission. The selection of self-gravitating cores from the synthetic continuum map is 
discussed in § 4 and the observational results are summarized in § 5. In § 6 the Chandrasekhar-Fermi 
formula is tested on the selected cores. Results are discussed in § 7, and conclusions are summarized in § 8. 



2. The MHD Model 

We solve the compressible MHD equations: 

<91np 



dt 



+ v Vlnp= -V-v, (1) 



dv PI 

— + V -Vv = VlnP+-j xB-V$ + f, (2) 

at p p 



P = pT, (3) 



3D 

— = V x (v x B), (4) 



j = VxB, 



(5) 
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V 2 $ = C*p 



(6) 



plus numerical diffusion terms, and with periodic boundary conditions, v is the velocity, B the magnetic 
field, $ the gravitational potential, f an external random force, T the gas temperature and T = const 
is assumed. The isothermal approximation is appropriate for the short cooling time in relatively dense 
molecular gas (see discussion in Padoan, Zweibel & Nordlund 2000). The constant C is given by: 



where the velocity is measured in units of vq, length in units of Iq, time in units of Iq/vq, and density in 
units of po- 

For the purpose of this paper we have run a numerical simulation of super-sonic, super-Alfvenic and 
self-gravitating MHD turbulence, by solving numerically the equations above, in a 128 3 mesh. The initial 
density and magnetic fields are uniform; the initial velocity is random, generated in Fourier space with 
power only on the large scale. We also apply an external random force, to drive the turbulence at a roughly 
constant rms Mach number of the flow. This force is generated in Fourier space, with power only on small 
wave numbers (1 < k < 2), as the initial velocity. 

We let the flow evolve for one dynamical time, defined as td yn — Lo/a v , where L is the linear size of 
the computational box, and a v is the rms flow velocity. A self-gravitating flow is never statistically relaxed. 
First a number of collapsing and accreting cores is generated, and later the cores begin to merge. The 
numerical resolution (128 3 numerical mesh) allows only the description of the initial phase of the collapse 
of single cores. Results are progressively inaccurate at later times, when the numerical resolution cannot 
cope with the exceedingly high density. We therefore interrupt the simulation at a time when most cores 
are just recently formed and start to collapse, which is about one dynamical time of the large scale. 

Periodic boundary conditions and large scale external forcing are justified by the fact that we simulate 
a region of turbulent flow inside a larger turbulent molecular cloud. The rms Mach number of the flow is 
M s ~ 11, which corresponds to a linear size L w 6.25 pc, and an average gas density < n >« 320 cm -3 , 
using empirical Larson type relations (Larson 1981). The average magnetic field in this model is rather 
weak, as justified by our previous work (Padoan & Nordlund 1997, 1999), and such that the average 
magnetic energy is approximately twice the thermal energy. Assuming a kinetic temperature T = 10 K, 
the rms flow velocity is a v = 3.0 km/s; the average magnetic field strength is (B) = 3.3 /iG, while 
(B 2 ) 1 / 2 = 7.1 fiG . Despite this low value of (B), strongly magnetized cores are formed by the process of 
turbulent fragmentation (the maximum field strength in the simulation is 127.8 /iG), due to compressions 
and stretching in the turbulent flow. 



We compute the polarized thermal dust emission from the MHD model following the formalism in Fiege 
& Pudritz (2000) (see also Wardle & Konigl 1990). As in Fiege & Pudritz (2000) we arc interested in the 
thermal dust emission at sub-mm wavelengths, and therefore the effect of self-absorption (Hildebrand et al. 
1999) and scattering (Novak et al. 1989) can be neglected. We further assume that the grain properties are 



C = 




(7) 
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Polarized Thermal Dust Emission 



- 5 - 



1.5 



Astronomical Silicate, X= 350yiim 








2 3 4 

axial ratio = (intermediate + long)/(3xshort) 



5 



Fig. 1. — Grain cross section ratio as a function of grain axial ratio for prolate spheroids (dashed line) and 
oblate spheroids (solid line), at wavelength A = 350 /im. The grains are the "astronomical silicate" discussed 
by Draine & Lee (1984), as modified by Li & Draine (2001). 

constant^] and the temperature is uniform. With these assumptions, the Stokes parameters Q and U are 
proportional to the following integrals along the line of sight: 



where p is the gas density, ip the angle between the projection of B on the plane of the sky and the North, 
7 the angle between the local B vector and the plane of the sky. The polarization angle \ is then given by: 




(8) 




(9) 



tan 2\ = — ■ 

q 



(10) 



The degree of polarization P, defined as 



P 



VQ 2 + u 2 



(11) 



where I is the intensity, is given by: 



P = a 



yjq 2 + u 2 



(12) 



2 Fiege and Pudritz (2000) consider explicitly the contributions to the polarized emission from several grain species, which 
is not required here since we do not specify any particular set of grain properties. 
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with 

S = J pds, (13) 

and 

£2 = \ fp (cos 2 7-2/3) da, (14) 



The coefficient a in equation (12) is defined as 



a = RF^2L, (15) 

where R, F, C po i and C ran are defined by Lee & Draine (1985). C po \ is the grain polarization cross section, 
which differs for oblate spheroids, 

Cpoi — C±_ - C||, (16) 

and prolate spheroids. 

C pol = (C {] - C±)/2, (17) 

where Cm and Cj_ are the polarization cross sections for electric field parallel and perpendicular to the 
symmetry axis of the grain. C ran is the average cross section for randomly oriented grains, 

C ran = (2C ± + C||)/3. (18) 

R is the Rayleigh polarization reduction factor due to imperfect grain alignment, and F the polarization 
reduction factor due to the turbulent component of the magnetic field. In our case F = 1 because the 
three dimensional magnetic field is provided by the MHD model (assuming that no small scale structure 
is unresolved in the numerical solution). However, in this work we do not attempt to use specific models 
for the grain cross sections and for the reduction factors, since our aim is to investigate the effect of the 
magnetic field topology. We therefore adopt a constant value of a = 0.15, which is equivalent to assuming a 
reasonable value for the maximum degree of polarization, P max = 15%, since 

Praax = ": TZ ■ (19) 

1 — a/6 



In Figure 1 we have plotted C po i/C ran of astronomical silicate spheroids, for "perfect spinning 
alignment", as a function of the "axial ratio", defined as (a + b)/(2c), where a, b and c are the long, 
intermediate and short axis respectively. For prolate spheroids b = c; for oblate spheroids a — b. With this 
definition, the ratio C po i/C ra n is nearly the same, at given "axial ratio", for prolate and oblate spheroids. 
The astronomical silicate is that discussed by Draine & Lee (1984), as modified by Li & Draine (2001). 
For wavelengths A in the range 100 /im < A < 1mm, C po i/C ran is almost independent of A. Figure 1 is 
computed for A = 350 /im. The adopted value of a — 0.15 would require only moderate deviations from 
sphericity, provided the degree of alignment R is reasonably high. If for example R ~ 0.33, a reasonable 
value of approximately 1.3 for the axial ratio gives the assumed value of a. 

We have computed polarization maps from the MHD model in two cases. In the first case, or model 
A, it is assumed that grains are aligned independently of their visual extinction, Ay] in the second case, 
or model B, it is assumed that grains are not aligned at visual extinction larger than Ay t a — 3 mag. The 
intensity of the radiation field was computed with a Monte Carlo method. Photon packages were sent into 
the cloud from the background and the scatterings and absorption processes were simulated. In each cell 
the number of incoming photons was registered and the intensity relative to the background was used to 



Fig. 2. — Polarization map from the MHD model (model A). The length of the polarization vectors is 
proportional to the degree of polarization, with the longest vector corresponding to P = 14 %. Only one 
polarization vector every three computational cells is plotted. The contour map shows the sub-mm dust 
continuum intensity /. The position of the self-gravitating cores CI, C2 and C3 (see text in §4) is also 
shown. 

compute the effective values of Ay. The dust opacity was calculated from the relation between -/V(H 2 ) and 
Ay as given by Bohlin, Savage & Drake (1978). The albedo of the grains was 0.5 and the asymmetry factor 
of the scattering g=0.6. 

The polarization maps are plotted in Figure 2 for both model A and in Figure 3 for model B. The 
average magnetic field (volume-average over the whole computational box) is oriented in the vertical 
direction. The length of the polarization vectors is proportional to the degree of polarization, with the 
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Fig. 3. — As in Figure 1, but assuming unaligned grains above Av,o = 3 mag (model B). 

longest vector corresponding to P = 14 %. The contour map shows the sub-mm dust continuum intensity 
I. Model B shows systematically low values of the degree of polarization P in the brightest regions, because 
at Ay > 3 mag the dust contributes to the total flux, but not to the polarized flux (P is defined as the ratio 
of polarized and total flux). This does not happen in model A, with the exception of one region (core C3). 

There is no obvious correlation between P and I over the entire map, as shown in the scatter plots of 
Figure 4. P spans a range of values 0.3 %< P < 14 %, with (P) = 5.6 %. Such relatively large values of P, 
comparable to typically observed values, are usually found in models of super-Alfvenic turbulence. This is 
due to the fact that each line of sight is normally dominated by one particular dense filament (or edge-on 
sheet), where the magnetic field is rather uniform because it is amplified in the directions perpendicular to 
the shock compression that generates that filament (or sheet) . The power spectrum of turbulence contains 
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most of the energy on the largest scales, which constrains the magnetic field tangling, in the sense that the 
number of correlation lengths of the magnetic field can never be ^ 1 (see the discussion in § 7) . 




Fig. 4. — Left panel: Scatter plot of the degree of polarization versus the emission intensity from the entire 
map in Figure 2. Only the polarization vectors shown in Figure 2 are used, that is one position every three 
computational cells. The average degree of polarization is (P) = 5.6 %. Right panel: As left panel, but from 
the entire map in Figure 3, that is with unaligned grains above Ay.o — 3 mag. 



4. Self Gravitating Cores 

For our polarization study, we select the brightest cores in the intensity map shown in Figure 2, in 
a similar way as protostellar cores would be selected from a sub-mm dust continuum map. Since we are 
interested in relating the observable properties of the cores with their intrinsic three dimensional structure, 
we have limited our selection to intensity maxima that are also well defined in the three dimensional MHD 
data-cube. We have excluded the remaining two or three bright cores that are not easily identified in the 
original three dimensional data-cube, because of the contribution of multiple structures along the line of 
sight. 

We have selected three cores that we call CI, C2 and C3. Their position is shown on the maps 
in Figures 2 and 3. We have verified, using the MHD data-cube, that the cores CI, C2 and C3 are 
gravitationally bound (gravitational energy in excess of the kinetic energy of turbulence). While the core 
CI is still being assembled by a turbulent shock (this is also indicated by its position as a density peak in 
a strongly curved segment of a filament), the core C3 is collapsing with a surrounding flow characterized 
by accretion along numerous filaments converging towards the core. The magnetic field around the core is 
roughly aligned with these accretion filaments. In the core CI the magnetic field is roughly aligned with 
the filament where the core is located (it is amplified in that direction by the shock compression), and is 
therefore strongly curved inside the core. The core C2 seems to be an intermediate case between CI and 
C3. 

The collapsing/accreting configuration of the core C3 can be explained as follows. Once a core is 
formed by a shock as a density peak along a filament (or on a sheet), part of the gas is moving along the 
filament (or on the sheet) away from the core, because of the pressure gradient. However, as long as the 
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core becomes massive enough, its gravitational force becomes larger than the pressure gradient and the flow 
along the filament is reversed back into the core. In the simple case of a single shock the collapsing core will 
be accreting along two filaments (the two sides of the original filament), while in the case of the brightest 
cores, located at the intersection of a number of shocks, many accreting filaments can be found. 




Fig. 5. — Upper panels: Polarization maps as in Figure 2, but limited to the regions around the cores CI, 
C2 and C3 (see text). The polarization vectors are here plotted for one position every two computational 
cells, and the largest polarization vectors (found in the core CI) corresponds to P rj 11.5 %. Lower panels: 
As upper panels, but from the map in Figure 3, that is with unaligned grains above Ay^o = 3 mag. 

The region of the polarization map around each core is shown in the upper panels of Figure 5 for 
model A. The polarization vectors are here plotted for one position every two computational cells, and the 
largest polarization vectors (found in the core CI) corresponds to P ~ 11.5 %. The respective scatter plots 
of P versus I (also for model A) are shown in the left panels of Figure 6, where the intensity / is plotted 
in linear units, for a more direct comparison with published observational plots. The difference in the P—I 
plot between the cores CI and C3 is remarkable in model A. The plot for the core CI is characterized by a 
lower envelope with P increasing with increasing /, while the plot for the core C3 shows a very clear upper 
envelope with P decreasing with increasing I. The excursion of the values of P on the envelopes is from 
approximately 2% to 9% in the first case and from approximately 8% to 2% in the second case. The plot 
for the core C2 looks like an intermediate case. 

Should we conclude that the P—I plot is a diagnostic for the time evolution of a core, from an early 
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Fig. 6. — Left panels (model A): Scatter plots of the degree of polarization versus the emission intensity 
from the maps in Figure 5 (upper panels). Only the polarization vectors shown in Figure 5 are used, that is 
one position every two computational cells. Right panels (model B): same as left panels, but from the maps 
in the lower panels of Figure 5, that is assuming unaligned grains above Ay,o = 3 mag. 
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Fig. 7. — Left panel (model A): Scatter plot of the degree of polarization versus the emission intensity for 
the core CI, seen from three orthogonal points of view. Right panels (model B): same as left panels, but 
with unaligned grains above Av,o — 3 mag. 
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shock phase to a collapse phase? The integration of the Stokes parameters along different lines of sight 
shows that the P—I plot of the core C3 depends on the direction of the line of sight. Along certain directions 
P and / are virtually uncorrelated; along others P decreases with increasing / because the magnetic field is 
roughly parallel to the line of sight towards the center of the core, or because of a cancellation effect due 
to the direction of polarization within material in front of the core. The P—I plot is therefore sensitive to 
the viewing angle, and is not a diagnostic for the evolutionary state of the cores (at least with the present 
numerical resolution that allows only the description of the very early phase of the core collapse) . 

The situation is different when grains are assumed to be aligned only at visual extinction below 
Av.o = 3 mag, as in model B. The polarization maps of the cores in the case of model B are shown in 
the lower panels of Figure 5, and the corresponding P—I plots in the right panels of Figure 6. The upper 
envelope in the P—I plot, with P decreasing with increasing /, is always found, since at visual extinction 
above Ayo = 3 mag the dust contributes to the continuum intensity /, but not to the polarized flux (P is 
defined as the ratio of polarized and total flux -eq. (|TT|)). The direction of polarization is also affected by 
the lack of alignement of grains at large visual extinction. As an example, the region at the south of the 
core C2 shows a magnetic field direction roughly parallel to the elongation of the core in model A, while the 
magnetic field direction is almost perpendicular to the core elongation in model B. 

Figure 7 shows the P—I plots for the core CI seen from three different orthogonal directions, in model 
A (left panels) and B (right panels). These plots show that P always decreases with increasing / in model 
B, independent of the viewing angle. 

The histograms of the polarization angle \ m the maps of Figure 5 are plotted in Figure 8. For 
computing these histograms we have used all the polarization vectors available, that is one per computational 
cell. Since the maps are made of 20 x 20 pixels, we have used 400 polarization vectors for each histogram, 
while in Figures 4 and 5 only 10 x 10 = 100 vectors were used for clarity. The plots in Figure 7 show 
very broad distributions of \- The values of \ are limited to the interval [— 7r/2, 7r/2], but the reference 
direction, or the origin of the interval, can be shifted arbitrarily by an amount between and w. We have 
done so in order to minimize the dispersion of \i which is obtained by shifting the origin of the angle to the 
approximate location of the minimum in the histogram. In this way the dispersion o x has a robust physical 
meaning, indipendent of the particular local orientation of the average magnetic field relative to the North. 
No Gaussian fit of the histogram is attempted, and o~ x is simply defined as 

^H(x-(x» 2 > 1/2 (20) 

with the values of % defined in the interval that minimizes this expression, as explained above. Differences 
between models A and B are rather small. Results are listed Table 1. For both models A and B we find 
the average value of the rms \ to be a x — 0.60 ±0.11. 

5. Observational P—I Plots of "Quiescent" Cores 

In this work we define as "quiescent" any starless or low mass star forming core. We do not try to 
model regions of massive star formation, where the enhanced radiation field must have a considerable effect 
on the grain alignment (Draine & Weingartner 1996, 1997). 

A number of recent sub-mm polarization maps of quiescent cores can be used to study the correlation 
of the degree of polarization P and the sub-mm dust continuum intensity /. Davis et al. (2000) observed 
the Serpens low mass star forming region, which contains eight cores. Polarization maps of the three cores 
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Fig. 8. — Left panels (model A): Histograms of the polarization angle x from the polarization maps shown in 
Figure 5 (upper panels). All the available polarization vectors have been used, that is one per computational 
cell. Right panels (model B): same as left panels, but from the polarization maps in the lower panels of 
Figure 5, that is with unaligned grains above Ay,o = 3 mag. 
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Fig. 9. — Spectral maps of the J=1~0 transition of CS. The maps are located around the cores CI, C2 and 
C3, from top to bottom respectively. The regions shown in the maps correspond to the regions included in 
the maps of Figure 5. The velocity axis ranges from -2.76 to 2.76 km/s, and the antenna temperature axis 
from to 3 K. 

L1544, L183 and L43 were obtained by Ward- Thompson ct al. (2000) and three more cores, CB 26, CB 24 
and DC 253-1.6 (CG 30) were mapped by Henning et al. (2001). All these regions were observed at 850 /im, 
with the SCUBA bolometer array at the James Clerk Maxwell Telescope. In all of the 14 cores the value 
of P is found to decrease rather steeply with increasing /. Polarization vectors from the lower intensity 
regions around the cores are always much larger than the vectors close to the center of the cores. 

P-I plots are presented both in Davis et al. (2000) and in Henning et al. (2001). They are very 
similar to the P-I plots from our model B, presented in the previous section. The upper envelope of the 
plots is defined by values of P from 1-2% at the highest /, to 8-10% at the lowest values of /. The drop in 
P is roughly consistent with a constant polarized flux, where the polarized flux is defined as P x /. This 
suggests that a significant fraction of the dust emission is from grains that do not contribute significantly 
to the polarization, presumably unaligned grains. 

Model A cannot reproduce the same result for all cores, but only in some cases due to the fortuitous 
orientation of the magnetic field or to a cancellation effect also dependent on the orientation of the core 
relative to the line of sight. It is possible that a numerical simulation with higher resolution, able to follow 
the collapse of cores to larger densities, may generate a stronger cancellation effect. The gravitational 
collapse is expected to drag the lines of force of the magnetic field towards the core center, which could 
generate a cancellation effect. However, it is not clear how important this effect may be, once diluted over 
the size of the telescope beam. It is also difficult to imagine it should be independent of orientation, since 
the average magnetic field before the collapse defines a specific direction. 



6. Distribution of Polarization Angle and Magnetic Field Strength 

In linear transverse oscillations of a line of force (for example Alfven waves in a uniform medium) the 
ratio of the lateral and transverse velocities is equal to the deviation of the line of force from a straight line. 
This deviation of the field from a straight line is measured by the angle x between the straight line and the 
local line of force. If the transverse velocity is the Alfven velocity v a , and the lateral velocity is the rms 
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turbulcnt velocity a v , then 




(21) 



where M a is the Alfvenic Mach number. This is strictly valid only in linear theory. We can dehne a 
coefficient / that quantifies the deviation from the prediction of the linear theory: 



where the numerator is the rms Alfven velocity, averaged inside the volume under study. Equation (|2lj), or 
cquivalcntly the value / = 1, was first used to derive the magnetic field strength in the Galaxy by Davis 
(1951) and Chandrasekhar and Fermi (1953). This method was further discussed by Zweibel (1990) and 
Myers & Goodman (1991). It has also been tested on the results of MHD simulations of sub-Alfvenic 
turbulence by Ostriker, Stone & Gammie (2001). They find that for small dispersion of the polarization 
angle (a x < 0.44 radians) the Chandrasekhar-Fermi formula provides a good estimate of the plane of the 
sky field strength, if a coefficient / « 0.5 is used. Similar results are found by Heitsch et al. (2001), using 
numerical simulations of mildly super- Alfvenic turbulence. 

In protostellar cores and Bok globules the dispersion of the polarization angle is sometimes very large, 
and the Chandrasekhar-Fermi formula is not expected to hold. However, the formula is sometimes applied 
irrespective of the large value of a x (eg Itoh et al. 1999; Davis et al. 2000; Henning et al. 2001). This is 
perhaps justified by the expectation that the degree of topological complexity of the magnetic field, and so 
the value of a x , should be related in some way to the rms Alfvenic Mach number M. a . A full numerical 
study of the dependence of a x on M. a can be pursued with the MHD turbulence simulations, and will be 
presented elsewhere. For the purpose of this work we test the Chandrasekhar-Fermi formula on the cores 
CI, C2 and C3, and provide the values of the corresponding correction factor /. 

The turbulent velocity dispersion in protostellar cores and Bok globules is often estimated using the 
line width of the CS molecule. It is possible that the bulk of the CS emission come from an envelope around 
the core, and not from the core interior, where most of the dust emission is generated, at least for the 
densest cores. This complicates the application of the Chandrasekhar-Fermi formula and its interpretation. 
For this reason, we have tested the formula by defining the turbulent velocity as given by the CS spectra. 

Maps of synthetic spectra of the J=l-0 transition of CS are computed using a non-LTE Monte Carlo 
radiative transfer code (Juvela 1997, 1998), from the three dimensional density and velocity fields generated 
in the numerical MHD experiment. The method of computing synthetic spectra was presented in Padoan et 
al. (1998). We assume a linear size of the computational box L = 6.25 pc, a uniform temperature T — 10 K, 
an rms turbulent velocity <j v — 3.0 km/s and an average gas density (n) = 320 cm~ 3 . These values are used 
to scale the numerical variables into physical units. In the radiative transfer computation the full three 
dimensional density and velocity fields are used. The radiative transfer is a non-local problem, and so it is 
solved for the whole three dimensional computational box. 

Spectral maps are obtained for the same lines of sight that define the polarization maps. The CS 
spectral maps that match the polarization maps of Figure 5 are plotted in Figure 9. All the spectra are 
used to compute the rms velocity a v . This is defined as the rms velocity averaged along the spectral profile, 
without any Gaussian fit: 



where T(v) is the antenna temperature at the velocity channel v, and the result is averaged over the whole 




(22) 




(23) 
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0.73 


1.32 


a x [rad] 


0.48 


0.73 


0.37 
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0.71 


0.55 


0.73 


0.65 




(0.49) 


(0.70) 


(0.37) 


(0.64) 


(0.62) 


(0.71) 


(0.52) 


(0.74) 


(0.65) 


/ 


0.32 


0.48 


0.29 


0.43 


0.29 


0.46 


0.35 


0.63 


0.31 




(0.32) 


(0.46) 


(0.29) 


(0.46) 


(0.31) 


(0.46) 


(0.33) 


(0.64) 


(0.31) 



Table 1: Estimated parameters for the three cores CI, C2 and C3. From top to bottom: rms turbulent 
velocity from CS spectral maps; rms polarization angle; correction factor for the Chandrasekhar-Fermi 
formula (see text). Values in parenthesis are for the case of unaligned grains above Av,o = 3 mag. 



spectral map of each core. The average velocity (v) is defined as 

JvT(v) dv 

{V) - jT(v)dv (24) 

The values of a v for three orthogonal lines of sight are listed in Table 1. The average value is 
g v = 1.01 ± 0.19 km/s. 

In order to determine the value of the coefficient /, we finally need to compute the rms Alfven velocity 
in the cores. This is done directly in the three dimensional MHD data cube. Using the values of a x 
computed in the previous section, and the values of a v and rms v a computed here, we can now evaluate the 
coefficient /. Results are listed in Table 1. The average value of / is found to be / = 0.40 ± 0.11 for both 
models A and B. 



7. Discussion 

In order to describe the formation of protostellar cores we have simulated the dynamics of molecular 
cloud turbulence on a relatively large scale (« 6 pc). Simulating a large scale is necessary because the 
boundary and initial conditions of a numerical simulation (and much more of an analytic model) are always 
idealized in some way. For example, a numerical simulation or an analytic model describing the evolution 
and collapse of a single protostellar core can hardly be considered as models for the core formation, since 
the core is obviously formed by the ad hoc setup of the initial conditions. On the other hand, if a larger 
scale is simulated, the process of core formation is controlled by the gas dynamics on the large scale, and 
may even be rather insensitive to the specific initial and boundary conditions at the large scale. 

In the simulation used for this work the typical size of cores is less than one tenth of the size of 
the computational box, and so we argue that our periodic boundary conditions do not affect too much 
the formation and evolution of dense cores. The setup of the initial conditions is also unlikely to affect 
dramatically the outcome of the simulation, since after one dynamical time the turbulent flow begins to 
lose memory of its initial conditions. The random external driving force can instead have important effects, 
in the sense that the flow would develop much more small scale structure if the forcing was done on a small 
scale, and the number of correlation lengths of the magnetic field, N corri could be relatively large. However, 
that would imply a power spectrum of turbulence peaking on small scales. In contrast, all observational 
evidence for the power spectrum of turbulence in molecular clouds is consistent with a scale free dynamics 
(Miesch & Bally 1994; Ossenkopf & Mac Low 2000), with most energy on the largest scale. The large scale 
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driving used in the present simulation, is therefore consistent with the observations. As a consequence, the 
number of correlation lengths of the magnetic field must be relatively small, N corr ~ 1. This explains the 
relatively large values of P found from our super-Alfvenic simulations, since a significant cancellation effect 
requires both a weak magnetic field (super-Alfvenic turbulence) and a large number of correlation lengths. 

It is well established that the polarized extinction of background stars at optical and near-IR 
wavelengths does not reveal the magnetic field direction in the dense interiors of molecular clouds, beyond 
an optical depth A v = 1 — 2 mag (Goodman et al. 1995; Arcc ct al. 1998). This is generally interpreted 
as the effect of unaligned grains at large visual extinction, since it can be shown that most alignment 
mechanisms are ineffective at large optical depth (Lazarian, Goodman & Myers 1997). Wiebe & Watson 
(2001), on the contrary, suggest that the results of polarization studies of background stars can be 
explained by a cancellation effect due to the magnetic field tangling along the line of sight. This requires 
super-Alfvenic turbulence and a large number of correlation lengths of the magnetic field, N corr . They find 
that N corr > 10 is necessary to explain the observations, and use a numerical MHD simulation by Stone, 
Ostriker & Gammie (1998) as an example of a theoretical model with N corr > 10. However, as explained in 
the previous paragraph, there is no observational evidence in favor of a power spectrum of turbulence with a 
peak at small scale, which is needed to produce a large value of N corr . The simulation by Stone, Ostriker & 
Gammie (1998) used by Wiebe and Watson (2001) is performed with a random driving on an intermediate 
scale, which explains the large value of N corr , but is certainly not motivated by the observations. 

In § 5 we have seen that recent sub-mm polarization maps of quiescent cores always produce P—I 
plots characterized by an upper envelope with P decreasing with increasing /. This result is predicted 
by our model B, that is by assuming that grains at visual extinction larger than Ay.o = 3 mag arc not 
aligned. It is not reproduced by model A, where the grain alignment is assumed to be independent of Ay, 
because the magnetic field tangling in the cores is insufficient to produce such drop in P. Although it is still 
possible that higher resolution simulations following the core collapse to higher densities generate a larger 
cancellation effect, that remains to be proven. From the present calculations it seems that the reasonable 
assumptions of model B provide a good interpretation of the observational P-I plots. If this interpretation 
is correct, we can conclude that sub-mm polarization maps cannot probe the magnetic field structure inside 
quiescent cores, at visual extinction larger than Ay = 3 mag. The use of such maps to constrain models 
of protostellar core formation and evolution is therefore questionable, especially if the effect of unaligned 
grains is not taken carefully into account. Polarization studies of sub-mm dust continuum emission from 
quiescent cores are therefore not inconsistent with previous results of polarized extinction of background 
stars at optical and near-IR wavelengths mentioned above. 

The specific value of Ay = 3 mag used to compute model B is a "conservative" choice. A value of 
Avfi ~ 2 mag, or even lower, is not ruled out by our study. It is likely that the alignment efficiency of 
grains begins to decrease significantly at Ay = 1 — 2 mag, as suggested by theoretical studies (Lazarian, 
Goodman & Myers 1997). In this work we show that at values of visual extinction larger then 3 mag grains 
must be almost completely unaligned. 

The MHD simulation used in this work has a limited numerical resolution. A 128 3 numerical mesh has 
been used and the volume selected around each core contains only 20 3 computational cells. Differences in 
the polarization maps of the cores between models A and B (see Figure 5) are equally limited. A much 
higher numerical resolution could in principle generates even larger differences between the magnetic field 
orientation in the densest regions of the cores and in their lower density envelopes. The center of the core 
at Ay > 3 mag, for example, could have collapsed dragging in the magnetic field lines, and this process 
would be totally invisible in the polarization maps. However, models for the formation and evolution of 



- 19 - 



protostellar cores could be tested with polarization maps of quiescent cores if the comparison were limited 
to the regions of low visual extinction around the cores. 

We have shown in the previous section that the Chandrasekhar-Fermi formula predicts the value of 
the rms Alfven speed in the cores, if a correction factor / « 0.4 is used, and if the rms turbulent velocity is 
derived from the line width of the J=l CS transition. The value of a x should be appropriately minimized 
as discussed in the previous section. However, in order to infer a magnetic field strength from an rms Alfven 
speed and a x , the gas density in the cores must be estimated, which introduces a significant uncertainty. 

8. Conclusions 

In this work we have computed dust continuum polarization from protostellar cores, assembled by 
super-sonic turbulent flows in models of molecular clouds, using numerical simulations of highly super-sonic 
MHD turbulence. The results are compared with recent polarization maps of quiescent protostellar cores 
and Bok globules. 

The main results of this work are: i) Values of P between 1 and 10% (up to almost P ma x) are produced 
by super-Alfvenic turbulence; ii) A steep decrease of P with increasing /, as observed in a number of 
protostellar cores and Bok globules, is always found in self gravitating cores selected from the MHD 
simulation, if grains are not aligned above a certain value of visual extinction Ay.o (model B); iii) The 
same behavior is hard to reproduce if dust grains are aligned independently of Ay (model A); iv) The 
Chandrasekhar-Fermi formula, corrected by a factor / « 0.4, can be used to estimate the average magnetic 
field strength in the cores. 

We conclude that sub-mm polarization maps of quiescent cores are not likely to faithfully map the 
magnetic field inside the cores at visual extinction larger than Ay.o ~ 3 mag. The use of such maps to 
constrain models of protostellar core formation and evolution is questionable, especially if the effect of 
unaligned grains is not properly taken into account. This conclusion shows that there is no inconsistency 
between the results from optical and near-IR polarized absorption of background stars, and the observed 
polarization of sub-mm dust continuum from quiescent cores. In both cases, grains at large visual extinction 
appear to be virtually unaligned. 

This work was supported by NSF grant AST-9721455. Ake Nordlund acknowledges partial support 
by the Danish National Research Foundation through its establishment of the Theoretical Astrophysics 
Center. Mika Juvela acknowledges support by the Academy of Finland Grant no. 1011055. 

References 

Arce, H. G., Goodman, A. A., Bastion, P., Manset, N., Sumner, M. 1998, ApJ, 499, L93 
Bohlin, R. C, Savage, B. D., Drake, J. F. 1978, ApJ, 224, 132 
Chandrasekhar, S., Fermi, E. 1953, ApJ, 118, 113 

Clemens, D. P., Bookbinder, J., Goodman, A., Kristen, H., Myers, P., Padoan, P., Wood, K., Heyer, M. H., 
Heiles, C, Jones, T. J., Dickey, J., Young, E., Rieke, G., Dow, K., Dowell, C. D., Drainc, B., Greaves, 



-20- 



J., Klaas, U., Laureijs, R., Lazarian, A., Shulz, B., Zweibel, E. 2000, in American Astronomical Society 
Meeting, Vol. 196, 2508 

Coppin, K. E. K., Greaves, J. S., Jenness, T., Holland, W. S. 2000, A&A, 356, 1031 

Davis, C. J., Chrysostomou, A., Matthews, H. E., Jenness, T., Ray, T. P. 2000, ApJ, 530, L115 

Davis, L. 1951, Phys. Rev., 81, 890 

Draine, B. T., Lee, H. M. 1984, ApJ, 285, 89 

Draine, B. T., Weingartner, J. C. 1996, ApJ, 470, 551 

Draine, B. T., Weingartner, J. C. 1997, ApJ, 480, 633 

Fiege, J. D., Pudritz, R. E. 2000, ApJ, 544, 830 

Gerakines, P. A., Whittet, D. C. B., Lazarian, A. 1995, ApJ, 455, L171 

Glenn, J., Walker, C. K., Young, E. T. 1999, ApJ, 511, 812 

Goodman, A. A., Jones, T. J., Lada, E. A., Myers, P. C. 1995, ApJ, 448, 748 

Greaves, J. S., Holland, W. S., Minchin, N. R., Murray, A. G., Stevens, J. A. 1999, A&A, 344, 668 



Heitsch, F., Zweibel, E., Mac Low, M. M., Li, P. S., Norman, M. L. 2001, |astro-ph/0103286 
Henning, T. 2001, ApJ, submitted 

Hildebrand, R. H., Dotson, J. L., Dowell, C. D., Schleuning, D. A., Vaillancourt, J. E. 1999, ApJ, 516, 834 
Itoh, Y., Chrysostomou, A., Burton, M., Hough, J. H., Tamura, M. 1999, MNRAS, 304, 406 
Juvela, M. 1997, A& A, 322, 943 
Juvela, M. 1998, A&A, 338, 723 
Larson, R. B. 1981, MNRAS, 194, 809 

Lazarian, A., Goodman, A. A., Myers, P. C. 1997, ApJ, 490, 273 
Lee, H. M., Draine, B. T. 1985, ApJ, 290, 211 
Li, A., Draine, B. T. 2001, ApJ, 553, (in press) 
Matthews, B. C, Wilson, C. D. 2000, ApJ, 531, 868 
Miesch, M. S., Bally, J. 1994, ApJ, 429, 645 

Minchin, N. R., Bonifacio, V. H. R., Murray, A. G. 1996, A&A, 315, L5 
Myers, P. C, Goodman, A. A. 1991, ApJ, 373, 509 

Novak, G., Gonatas, D. P., Hildebrand, R. H., Piatt, S. R., Dragovan, M. 1989, ApJ, 345, 802 



Ossenkopf, V., Mac Low, M. 2000, |astro-ph/0012247 

Ostriker, E. C, Stone, J. M., Gammie, C. F. 2001, ApJ, 546, 980 

Padoan, P., Juvela, M., Bally, J., Nordlund, A. 1998, ApJ, 504, 300 



- 21 - 



Padoan, P., Juvela, M., Goodman, A. A., Nordlund, A. 2001, ApJ, (in press -also |astro-ph/001112^ ) 

Padoan, P., Nordlund, A. 1997, |astro-ph/970617£ 

Padoan, P., Nordlund, A. 1999, ApJ, 526, 279 

Padoan, P., Zweibel, E., Nordlund, A. 2000, ApJ, 540, 332 

Vallee, J. P., Bastien, P., Greaves, J. S. 2000, ApJ, 542, 352 

Ward-Thompson, D., Kirk, J. M., Crutcher, R. M., Greaves, J. S., Holland, W. S., Andre, P. 2000, ApJ, 
537, L135 

Wardle, M., Konigl, A. 1990, ApJ, 362, 120 
Wiebe, D. S., Watson, W. D. 2001, ApJ, 549, L115 
Zweibel, E. G. 1990, ApJ, 362, 545 



This preprint was prepared with the AAS IATj^X macros v4.0. 



